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CONVERGENCE RATE FOR A GAUSS COLLOCATION METHOD 
APPLIED TO UNCONSTRAINED OPTIMAL CONTROL * 

WILLIAM W. HAGERt, HONGYAN HOU*, AND ANIL V. RAO§ 

Abstract. A local convergence rate is established for an orthogonal collocation method based on 
Gauss quadrature applied to an unconstrained optimal control problem. If the continuous problem 
has a sufficiently smooth solution and the Hamiltonian satisfies a strong convexity condition, then 
the discrete problem possesses a local minimizer in a neighborhood of the continuous solution, and 
as the number of collocation points increases, the discrete solution convergences exponentially fast in 
the sup-norm to the continuous solution. This is the first convergence rate result for an orthogonal 
collocation method based on global polynomials applied to an optimal control problem. 

Key words. Gauss collocation method, convergence rate, optimal control, orthogonal colloca¬ 
tion 


1. Introduction. A convergence rate is established for an orthogonal collocation 
method applied to an unconstrained control problem of the form 

minimize C'(x(l)) 

subject to x(t) = f(x(t), u(t)), 1,1], (1.1) 

x(-l) = Xo, 

where the state x(t) G R", x = ^x, the control u(t) G R™, f : R" x R™ —>• R", 
C : R" —>■ R, and Xp is the initial condition, which we assume is given. Assuming the 
dynamics x(t) = f(x(t), u(t)) is nice enough, we can solve for the state x as a function 
of the control u, and the control problem reduces to an unconstrained minimization 
over u. 

Let Vjsi denote the space of polynomials of degree at most N defined on the 
interval [—1, +1], and let denote the n-fold Cartesian product x ... x Pjy. We 
analyze a discrete approximation to of the form 

minimize C'(x(l)) 

subject to x(ri) = f(x(Ti), Ui), l<i<N, (1.2) 

x(-l) = Xp, X G 

The collocation points Tj, 1 < i < N , are where the equation should be satisfied, and 
Ui is the control approximation at time t^. The dimension of Pat is A +1, while there 
are N + 1 equations in () 1 . 2 () corresponding to the collocated dynamics at N points 
and the initial condition. When the discrete dynamics is nice enough, we can solve 
for the discrete state x G as a function of the discrete controls u^, 1 < i < A, 
and the discrete approximation reduces to an unconstrained minimization over the 
discrete controls. 
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We analyze the method developed in B HU where the collocation points are 
the Gauss quadrature abscissas, or equivalently, the roots of a Legendre polynomial. 
Other sets of collocation points that have been studied include the Lobatto quadrature 
points 011], the Chebyshev quadrature points Bi, the Radau quadrature points 
[ini[iii[Ta[2T], and extrema of Jacobi polynomials [24] • The Gauss quadrature points 
that we analyze are symmetric about t = 0 and satisfy 

— 1<Ti<T2<...<Tjv< +1. 

In addition, we employ two noncollocated points tq = —1 and tjv+i = +1- 

Our goal is to show that if (x*, u*) is a local minimizer for then the discrete 

problem (11.21) has a local minimizer (x^,u^) that converges exponentially fast in 
N to (x*,u*) at the collocation points. This is the first convergence rate result for 
an orthogonal collocation method based on global polynomials applied to an optimal 
control problem. A consistency result for a scheme based on global polynomials and 
Lobatto collocation is given in m- Convergence rates have been obtained previously 
when the approximating space consists of piecewise polynomials as in lUlllsilKIl 
(IHIHU- In these results, convergence is achieved by letting the mesh spacing tend 
to zero. In our results, on the other hand, convergence is achieved by letting iV, the 
degree of the approximating polynomials, tend to infinity. 

To state onr convergence results in a precise way, we need to introduce a function 
space setting. Let C^(R”) denote the space of k times continuously differentiable 
functions x : [—1, +1] —)► R" with the sup-norm || • ||oo given by 

||x||oo = sup{|x(t)| : t G [-1, -pl]}, (1.3) 

where | • | is the Euclidean norm. It is assumed that (|l.ll) has a local minimizer (x*, u*) 
in C^(K") X Given y G K", the ball with center y and radius p is denoted 

Sp(y) = {xGK":|x-y|<p}. 

It is assumed that there exists an open set C and p > 0 such that 

Bp{:x*(t),u*(t)) C n for all t G [—1, -|-1]. 

Moreover, the first two derivative of / and C are continuous on the closure of fl and 
on Sp(x*(l)) respectively. 

Let A* denote the solntion of the linear costate equation 

A*(t) = -V,iJ(x*(t),u*(t),A*(t)), A*(l)=VC(x*(l)), (1.4) 

where H is the Hamiltonian defined by iJ(x, u. A) = A^f(x, u). Here VC denotes 
the gradient of C. By the first-order optimality conditions (Pontryagin’s minimum 
principle), we have 


V„iJ(x*(t),u*(t),A*(t)) =0 


(1.5) 


for all t G [—1, -|-1]. 

Since the discrete collocation problem (II.2p is finite dimensional, the first-order 
optimality conditions (Karush-Kuhn-Tucker conditions) imply that when a constraint 
qualification holds [50], the gradient of the Lagrangian vanishes. By the analysis in 
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m, the gradient of the Lagrangian vanishes if and only if there exists A £ such 
that 


A(ri) =(x(ri),Uj, A(ri)), 1 < i < iV, (1.6) 

A(+l) = VC(x(+l)), (1.7) 

0 = (x(Ti),Ui, A(tj)) , 1 < i < A^. (1.8) 

The assumptions that play a key role in the convergence analysis are the following: 

(Al) X* and A* £ for some rj > 3. 

(A2) For some a > 0, the smallest eigenvalue of the Hessian matrices 
V2C(x*(l)) and V^,_,)i7(x*(t),u*(t),A*(t)) 

is greater than a, uniformly for t £ [—1,+!]. 

(A3) The Jacobian of the dynamics satisfies 

||V,f(x*(t),u*(t))|U < 1/4 and || V,f(x*(t), u*(t))T||^ < 1/4 

for all t £ [—1, +1] where || • ||o<d is the matrix sup-norm (largest absolute row 
sum), and the Jacobian Va;f is an n by n matrix whose i-th row is ■ 

The smoothness assumption (Al) is used to obtain a bound for the accuracy with 
which the interpolant of the continuous state x* satisfies the discrete dynamics. The 
coercivity assumption (A2) ensures that the solution of the discrete problem is a local 
minimizer. The condition (A3) does not appear in convergence analysis for (local) 
piecewise polynomial techniques [3ii[5iia[ni[Ta[22]. It arises when we approximate 
a solution by polynomials defined on the entire interval [—1,-|-1]. More precisely, in 
the analysis, the dynamics is linearized around (x*,u*), and (A3) implies that when 
we perturb the linearized dynamics, the state perturbation is bounded uniformly 
in N with respect to the perturbation in the dynamics. If the domain [—1,-|-1] is 
partitioned into uniform subdomains of width h and a different polynomial is used on 
each subdomain, then (A3) is replaced by 

||V,f(x*(t),u*(t))|U < l/(2h) and ||V,f(x*(t),u*(t))T||^ < l/(2/i) 

which is satisfied when h is sufficiently small. In general, (A3) could be replaced by 
any condition that ensures stability of the linearized dynamics. 

In addition to the 3 assumptions, the analysis employs 2 properties of the Gauss 
collocation scheme. Let ujj, 1 < j < N, denote the Gauss quadrature weights, and 
for 1 < z < and 0 < j < N , define 


N 

Dij = Lj{n), where Lj{t) := . (1.9) 

■ n 

D is a differentiation matrix in the sense that (Dp)^ = piri), 1 < * < A^, where 
p £ Vn is the polynomial that satisfies p{Tj) = pj for 0 < j < N. The submatrix 
Di ■,N consisting of the tailing N columns of D, has the following properties: 

(PI) T)i n is invertible and ||Dj^.]y^||oo < 2. 

(P2) If W is the diagonal matrix containing the quadrature weights on the 
diagonal, then the rows of the matrix [W^/^Di: 7 v]~^ have Euclidean norm 
bounded by v^. 
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The fact that Di is invertible is established in [H Prop. 1]. The bounds on the 
norms in (PI) and (P2), however, are more subtle. We refer to (PI) and (P2) as 
properties rather than assumptions since the matrices are readily evaluated, and we 
can check numerically that (PI) and (P2) are always satisfied. In fact, numerically we 
find that ||D)f]y||oo = 1 Ttat where the last Gauss quadrature abscissa rjv approaches 
+1 as N tends to oo. On the other hand, we do not yet have a general proof for the 
properties (PI) and (P2). In contrast, conditions (AI)-(A3) are assumptions that are 
only satisfied by certain control problems. 

If G is a solution of (11.21) associated with the discrete controls Uj, 1 < z < 



satisfies (ll.6l)-(ll.8l). then we 

define 


x^ = 

[ x^(-l), 

x^(ti), ..., 

X^(rAr), 

x^(+l) 

X* = 

[ x*(-l), 

x*(ri), ..., 

X*(rAr), 

x*(+l) 

= 

[ 

Ui, 

Uat 


u* = 

[ 

u*(ti), ..., 



A^ = 

[ A^(-l), 

A"(ri), ..., 

A^(rAr), 

A^(+l) 

A* = 

[ A*(-l), 

A*(ti), ..., 

A* (tat), 

A*(+l) 


For any of the discrete variables, we define a discrete sup-norm analogous to the 
continuous sup-norm in (11.311 . For example, if G with Ui G R™, then 

||U'^||oo=sup{|U,|:l<z<iV}. 

The following convergence result is established: 

Theorem 1.1. If (x*, u* ) is a local minimizer for the continuous problem dnj 
and both (A1)-(A3) and (P1)-(P2) hold, then for N sufficiently large with N > rj + l, 
the discrete problem dill) has a local minimizer (X^,U^) and an associated discrete 
costate for which 

max {IIX^^ - X*||oo, IIU"^ - U*||oo, - AIU} < cN^-^ (1.10) 

where c is independent of N. 

Although the discrete problem only possesses discrete controls at the collocation 
points —1 < Ti < -1-1, 1 < z < A^, an estimate for the discrete control at t = 
— 1 and t = -|-1 is usually obtained from the minimum principle (11.51) since we do 
have estimates for the discrete state and costate at the end points. Alternatively, 
polynomial interpolation could be used to obtain estimates for the control at the end 
points of the interval. 

The paper is organized as follows. In Section [H the discrete optimization problem 
dEH) is reformulated as a nonlinear system of equations obtained from the first-order 
optimality conditions, and a general approach to convergence analysis is presented. 
Section [3] obtains an estimate for how closely the solution to the continuous prob¬ 
lem satisfies the first-order optimality conditions for the discrete problem. Section |4] 
proves that the linearization of the discrete control problem around a solution of the 
continuous problem is invertible. Section [5] establishes an stability property for 
the linearization, while Section |6] strengthens the norm to L°°. This stability prop¬ 
erty is the basis for the proof of Theorem 11.11 A numerical example illustrating the 
exponential convergence result is given in Section H 

Notation. The meaning of the norm || • ||oo is based on context. If x G C°(R”), 
then ||x||oo denotes the maximum of |x(t)| over t G [—1, -1-1], where | • | is the Euclidean 
norm. If A G R"*^", then ||A||oo is the largest absolute row sum (the matrix norm 
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induces by the £oo vector norm). If U € is the discrete control with XJi € K."*, 

then ||U||oo is the maximum of |Ui|, 1 < i < The dimension of the identity matrix 
I is often clear from context; when necessary, the dimension of I is specihed by a 
subscript. For example, I„ is the n by n identity matrix. VC denotes the gradient, 
a column vector, while V^C denotes the Hessian matrix. Throughout the paper, c 
denotes a generic constant which has different values in different equations. The value 
of this constant is always independent of N. 1 denotes a vector whose entries are all 
equal to one, while 0 is a vector whose entries are all equal to zero, their dimension 
should be clear from context. 

2. Abstract setting. As shown in [T^], the discrete problem can be refor¬ 
mulated as the nonlinear programming problem 

minimize C(Xjv+i) 

subject to = f(Xj,Uj), 1 < i < Xq = Xq, (2.1) 

X^+I = Xo + Ef=i (X„ U,). 

As indicated before Theorem ll.il X^ corresponds to x^(ri). Also, [12] shows that the 
equations obtained by setting the gradient of the Lagrangian to zero are equivalent 
to the system of equations 


Af-l-l 


E 

i^J^-A, =-V,i/(X„U„A,), 

l<i<N, 

An+i = VC{Xn+i), (2.2) 

i=i 

0 = VuH (Xi, Ui, Ai) , 

l<i<N, 

(2.3) 

where 


l<i<N, 

l<j<N, (2.4) 


N 




q 

WI 

1 

II 

+ 

<i< N. 

(2.5) 


i=i 


Here A^ corresponds to A'^(ri). The relationship between the discrete costate A^, the 
KKT multipliers associated with the discrete dynamics, and the multiplier Ajv+i 
associated with the equation for X^v+i is 

uiiAi = \i + u!i\N+i when 1 < i < A^, and Aat+i = Aat+i. (2.6) 


The first-order optimality conditions for the nonlinear program (EH) consist of 
the equations (12.21) and (12.31) . and the constraints in (12.11) . This system can be written 
as T(X, U, A) = 0 where 


(ri,r2,...,r5)(x,u,A) e x r" x r"^ x r" x r™^. 
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The 5 components of T are defined as follows: 

rH(x,u,A) = AjX.j - f(x„u,), 1 < ^ < iv, 

N 

r 2 (X, U, A) = X^+i - Xo - ^cc,f(X,, U,), 

i=i 

r 3 ,(X,U,A)= +V,77(X,,U„A,), 1 < ^ < A, 

r4(X, U, A) = A^+i - V,C(X^+i), 
r5.(X,U,A) = V„A(X„U„A,), 1 <Z<A. 

Note that in formulating T, we treat Xq as a constant whose value is the given 
starting condition xq. Alternatively, we could treat Xq as an unknown and then 
expand T to have a 6-th component Xq — Xq. With this expansion of T, we need 
to introduce an additional multiplier Aq for the constraint Xq — xq. To achieve a 
slight simplification in the analysis, we employ a 5-component T and treat Xq as a 
constant, not an unknown. 

The proof of Theorem 11.11 reduces to a study of solutions to T(X, U, A) = 0 in 
a neighborhood of (X*, U*, A*). Our analysis is based on 0 Proposition 3.1], which 
we simplify below to take into account the structure of our T. Other results like this 
are contained in Theorem 3.1 of [3], in Proposition 5.1 of [H, and in Theorem 2.1 of 

m 

Proposition 2.1. Let X be a Banaeh space and y be a linear normed space with 
the norms in both spaces denoted || ■ ||. Let T: X i —> y with T continuously Frechet 
dijferentiable in Br{d*) for some 6* £ X and r > 0. Suppose that 

||vr(0) - VT{e*)\\ < e for all 6 £ Br{e*) 

where VT(0*) is invertible, and define p := || VT(0*)“^||. If £p < 1 and ||T (0*)|| < 
(1 — p£)r/p, then there exists a unique 6 £ Br{d*) such that T{9) = 0. Moreover, we 
have the estimate 

\\e-e*\\<,^\\ri9*)\\<r. ( 2 . 7 ) 

1 — pe 

We apply Proposition 12.11 with 9* = (X*,U*,A*) and 9 = (X^,U^, A^). The 
key steps in the analysis are the estimation of the residual ||T(0*)||, the proof that 
VT{9*) is invertible, and the derivation of a bound for ||VT(0*)“^|| that is indepen¬ 
dent of N. In our context, for the norm in X, we take 

||0|| = ||(X,U, A)|U = max{||X|U, |1U|U, ||A|U}. (2.8) 

For this norm, the left side of (11.101) and the left side of (12.711 are the same. The 
norm on y enters into the estimation of both the residual ||T(0*)|| in (12.711 and the 
parameter p := ||VT(0*)“^||. In our context, we think of an element of 3^ as a vector 
with components y^, I < * < 3A -|- 2, where £ R" for I < * < 2 A -|- 2 and yi £ R™ 
for i > 2A -I- 2. For example, 7i(X,U, A) £ R"^ corresponds to the components 
Yi £ R", I < i < A. For the norm in 3^, we take 

||y||oo = sup{|yj| : I < z < 3A -b 2}. 


(2.9) 
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3. Analysis of the residual. We now establish a bound for the residual. 
Lemma 3.1. If (Al) holds, then there exits a constant c independent of N such 

that 


||r(X*,U*,A*)|U<c7V2-’’ 


(3.1) 


for all N > rj + 1. 

Proof. By the dehnition of T, 74 (X*,U*,A*) = 0 since x* and A* satisfy the 
boundary condition in (11.41) . Likewise, 75 (X*,U*, A*) = 0 since dni) holds for all 
t € [—1,+!]; in particular, (II. 5|) holds at the collocation points. 

Now let us consider 7i. By [TU Eq. (7)], 

N 

1=0 

where G is the (interpolating) polynomial that passes through x*(Ti) for 0 < 
i < N. Since x* satishes the dynamics of (II.IL it follows that f(X*,U*) = x*(Ti). 
Hence, we have 


rn(X*,U*,A*) =x^(r,)-x*(T,). 


(3.2) 


We combine Proposition 2.1 and Lemma 2.2 in m to obtain for A > 77 + 1 , 


llx^ - X* 


6 e 


A- 1 


(1 + 2A2) + 6 eA(l + ci^A) 


/ I2||a;(’^+^)| 
J 77 + 1 


where is the (77 + l)-st derivative of x and ci'/N is a bound for the Lebesgue 

constant of the point set n, 0 < i < N, given in Theorem 4.1 of [16]. Hence, there 
exists a constant C 2 , independent of A but dependent on 77 , such that 

||x^ - x*||oo < C 2 A^“''. (3.3) 


Consequently, 7i(X*,U*, A*) complies with the bound (13.11) . 
Next, let us consider 


N 

r2(X*,U*,A*)=x*(l)-x*(-l)-^c.,f(x*(r,),u*(r,)). (3.4) 

1=1 

By the fundamental theorem of calculus and the fact that A-point Gauss quadrature 
is exact for polynomials of degree up to 2A — 1, we have 

0 = {1) - {-!) - x^{t)dt = x^{l)-x^{-l)-J2^j^^i'^j)- (3-5) 

j=i 

Subtract (|3.5|) from (|3.4p to obtain 

N 

r 2 (X*, U*, A*) = X* (1) - x^ (1) + ^ a;, (x^ (r,) - x* (r,)) 

1=1 


(3.6) 
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Since uji > 0 and their sum is 2, it follows (13.31) that 

i=i 

By Theorem 15.1 in [53] and Lemma 2.2 and Theorem 4.1 in m, we have 
|X*(1)-X^(1)|<||X*-X^|U 

<(l+c.^jv)(^) (^)’’*‘||xC+‘>|U. (3.8) 

We combine (I3.6I) - (I3.8I) to see that 72(X*,U*, A*) complies with the bound (13.11) . 
Finally, let us consider Ti. By [121 Thm. 1], 

N + l 

^4A*=A^(rO, l<*<iV, 

j=i 

where G is the (interpolating) polynomial that passes through A* = A(tj) for 
1 < J < A + 1. Since A* satisfies (11.41) . it follows that A*(ri) = — V2;7L(X*, U*, A*). 
Hence, we have 

r3.(X*,U*,A*) = A^(r,)-A*(T,). 

Exactly as we handled 7i in (13.21) . we conclude that Ts^X*, U*, A*) complies with the 
bound (13.11) . This completes the proof. 0 

4. Invertibility. In this section, we show that the derivative VT{9*) is invert¬ 
ible. This is equivalent to showing that for each y G 3^, there is a unique 9 G X such 
that VT{9*)[9] = y. In our application, 9* = (X*,U*, A*) and 9 = (X,U,A). To 
simplify the notation, we let VT* [X, U, A] denote the derivative of T evaluated at 
(X*,U*, A*) operating on (X,U, A). This derivative involves the following 6 matri¬ 
ces: 

A,, = V,f(x*(rO,u*(TO), B, = V„f(x*(T0,W(T,)), 

Qi = '^xxH (x*(Ti),U*(Ti),A*(Ti)) , Si = Vuxff (x*(Ti),U*(Ti),A*(Ti)) , 

R, = V,„iL(x*(T,),u*(r,),A*(T,)), T = V2C(x*(l)). 


< 2c2JV^~^. 


(3.7) 


With this notation, the 5 components of VT*[X, U, A] are as follows: 


VTi* [X, U, A] = AjX, j - A,X, - B,U., 1 < i < JV, 

N 

vr 2 *[x, U, A] = Xiv+i - ^c^j(A,X, + B,U,), 

i=i 

VTg* [X, U, A] = ^ ^ j + AjA, + Q,X, -b S,U„ 1 < i < N, 

vr/ix, U, A] = Aiv+i - TXiv+ 1 , 

V%* [X, U, A] = Sjx, + R,U, + BjA,, l<t<N. 
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Notice that Xq does not appear in VT* since Xq is treated as a constant whose 
gradient vanishes. 

The analysis of invertibility starts with the first component of VT*. 

Lemma 4.1. If (PI) and (A3) hold, then for each q £ R” and p £ with 
Pi £ R", the linear system 

AjX, j - A,X, = p, l<i<N, (4.1) 

N 

X^+I - ^ c.,(A,X, + B,U,) = q, (4.2) 

has a unique solution X^- £ R”, 1 < J < A + 1. This solution has the bound 

||X,||oo<4||p|U + ||q||oo, l<j<N + l. (4.3) 


Proof Let X be the vector obtained by vertically stacking Xi through Xjv, let 
A be the block diagonal matrix with i-th diagonal block A^, 1 < i < N, and define 
D = Di:Ar 0 In where 0 is the Kronecker product. With this notation, the linear 
system (14.11) can be expressed (D —A)X = p. By (PI) Di:Ar is invertible which implies 
that D is invertible with Moreover, ||D“^||(x) = ||Di^]vllo° — ^ 

by (PI). By (A3) ||A|U < 1/4 and [iD-iAlU < ||D-i|U||A|U < 1/2. B_y [HI p. 
351], I — D~^A is invertible and ||(I — < 2. Consequently, D — A = 

D(I-D lA) is invertible, and 

||(D - A)-i||oo < 11(1 - D-iA)-i|U||D-i|U < 4. 

Thus there exists a unique X such that (D — A)X = p, and 


By gj), we have 


||X,|U<4||p|loo, l<j<N. 


(4.4) 


N 

||X^+i|U< ||q||oo + ^a;,||A,|U||X,|U. (4.5) 

i=i 

Since ||Aj||oo < 1/4 by (A3) and the ojj are positive and sum to 2, (14.41) and (14. 5|) 
complete the proof of (14.21) . □ 

Next, we establish the invertibility of VT*. 

Proposition 4.2. If (PI), (A2) and (A3) hold, then VT* is invertible. 

Proof. We formulate a strongly convex quadratic programming problem whose 
first-order optimality conditions reduce to VT*[X,U, A] = y. Due to the strong 
convexity of the objective function, the quadratic programming has a solution and 
there exists A such that VT* [X, U, A] = y. Since T* is square and VT* [X, U, A] = y 
has a solution for each choice of y, it follows that VT* is invertible. 

The quadratic program is 

minimize iQ(X, U) -L £(X, U) 

subject to AjXj = AiXi + BjUj -I- yii, l<i<N, ' 

Xn+ 1 = Ef=i (A,X, -f B,U,) + y2. 


(4.6) 
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where the quadratic and linear terms in the objective are 

N 

Q(X,U) = XT+iTXjv+i (XjQ.X, + 2 XTs,U, +UTr,U,) , (4.7) 

N 

C{X, U) = yJX ^+1 (yT X, + y^ U,) . (4.8) 

The linear term was chosen so that the first-order optimality conditions for (irei) 
reduce to VT*[X,U,A] = y. See [T2] for the manipulations needed to obtain the 
first-order optimality conditions in this form. By (A2), we have 


Q(X,U) >a 


N 


|X^+i|2 + ^u;, (|Xy+ |U, 




(4.9) 


Since a and lj are strictly positive, the objective of (14.61) is strongly convex, and by 
Lemma im the quadratic programming problem is feasible. Hence, there exists a 
unique solution to (14.61) for any choice of y, and since the constraints are linear, the 
first-order conditions hold. Consequently, VT*[X,U,A] = y has a solution for any 
choice of y and the proof is complete. □ 

5. w-norm bounds. In this section we obtain a bound for the (X, U) component 
of the solution to VT*[X, U, A] = y is terms of y. The bound we derive in this section 
is in terms of the w-norms defined by 


N 


N 


IXIl^ =|x^r+i|2 + ^a;.|xy and ||U||2 = ^ a;.|U,| 


(5.1) 


2=1 


2=1 


This defines a norm since the Gauss quadrature weight a;^ > 0 for each i. Since the 
(X, U) component of the solution to VT*[X, U, A] = y is a solution of the quadratic 
program (Irel) . we will bound the solution to the quadratic program. 

First, let us think more abstractly. Let tt be a symmetric, continuous bilinear 
functional defined on a Hilbert space "H, let £ be a continuous linear functional, let 
<() £ "H, and consider the quadratic program 

min I + (j),v + (f)) + £{v + (j>) : v GV 


where V is a subspace of 'H. If w is a minimizer, then by the first-order optimality 
conditions, we have 

n{w, v) + 7r(()), v) + i{v) = 0 for all v GV. 

Inserting v = w yields 


tt{w,w) =—{Tr{w, (p) + £{wj). (5-2) 

We apply this observation to the quadratic program (HU). We identify £ with 
the linear functional C in ()4.8p . and tt with the bilinear form associated with the 
quadratic term (14.7p . The subspace V is the null space of the linear operator in (14.6p 
and is a particular solution of the linear system. The complete solution of (14.61) is 
the particular solution plus the minimizer over the null space. 


GAUSS COLLOCATION 


11 


In more detail, let x denote the solution to (I4.1I) - (I4.2I) given by Lemma 14.11 for 
p = yi and q = y 2 - We consider the particular solution (X, U) of the linear system in 
(14.61) given by (x,0). The relation (15.21) describing the null space component (X,U) 
of the solution is 


Q(x,u) 


N 


{xn + y4)'^TXAr + V Wj [(Q*X* - yaO'^X, 


i=l 



(5.3) 


Here the terms containing x ^'I'e associated with TT{w,<f)), while the remaining terms 
are associated with i or equivalently, with C. By (A2) we have the lower bound 

Q(X,U)>a(||X|l2 +||U||2). (5.4) 


All the terms on the right side of (|5.3I) can be bounded with the Schwarz inequality; 
for example, 

N / N \ / Af 

^WiyJ.X, < I ^w.lyaip j I ^Wi|Xip 

i=l \i=l / \i=l 

<y2||y3||oo(||X|l2+||U|l2)'/^ (5.5) 

The last inequality exploits the fact that the oJi sum to 2 and |y 3 i| < Hyalloo- To 
handle the terms involving x in (15.31) . we utilize the upper bound ||Xj||oo < 5||y||oo 
based on Lemma 14.11 with p = yi and q = y 2 • Combining upper bounds of the form 
(j5.5() with the lower bound (15.4() . we conclude from (|5.3I) that both ||X|1^ and ||U||ij 
are bounded by a constant times ||y||oo- The complete solution of (14.6|) is the null 
space component that we just bounded plus the particular solution (x,0). Again, 
since ||Xj||oo < 5||y||oo, we obtain the following result. 

Lemma 5.1. If (A2)-(A3) and (PI) hold, then there exists a constant c, inde¬ 
pendent of N, such that the solution (X, U) of (14.61) satisfies ||X||ij < c||y||oo cind 
||U||^ < c||y||oo- 


6. oo-norm bounds. We now need to convert these w-norm bounds for X and 
U into oo-norm bounds and at the same time, obtain an oo-norm estimate for A. By 
Lemma mu the solution to the dynamics in (14.61) can be expressed 

X = (I-D”iA)-iD-iBU-Lp, (6.1) 


where B is the block diagonal matrix with i-th diagonal block B^. Taking norms and 
utilizing the bounds ||p||oo < 4||yi||oo and ||(I — D^^A)~^||^ < 2 from Lemma ITTI 
we obtain 


||X|U < 2||D-iBU||oo +4||yi||oo. (6.2) 

We now write 

D^^BU = [D-]v (g) I„]BU = [(Wi/2Di,Ar)-i 0 I„]BU^, (6.3) 

where W is the diagonal matrix with the quadrature weights on the diagonal and 
is the vector whose t-th element is y/uji iJi. Note that the y/uJi factors in (16.31) 
cancel each other. An element of the vector D^^BU is the dot product between 
a row of (W^/^Di: 7 v)“^ 0 In and the column vector BU,^. By (P2) the rows of 
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^ 0 I„ have Euclidean length bounded by v^. By the properties of 
matrix norms induced by vector norms, we have 

||BUJ|2< ||B||2||UJ|2 = ||B||2||U|U. 

It follows that 

IID-^BUIIoo < y2||B||2||U|U. (6.4) 

Combine Lemma EH] with (16.21) and (16.41) to deduce that ||X||oo < c||y||oo, where c 
is independent of N. Since |Xjv| < c||y||oo by Lemma ISHl it follows that ||X||oo < 
c||y||oo, 

Next, we use the third and fourth components of the linear system VT* [X, U, A] = 
y to obtain bounds for A. These equations can be written 

DtA + Djv+iAvr+i+ATA + QX + SU = y3 (6.5) 


and 


Aat+i — TXv+i = y4, (6.6) 

where A is obtained by vertically stacking Ai through Ajv, Q and S are block diagonal 
matrices with i-th diagonal blocks Qi and respectively, 0 I„, and 

D]y+i = 0 In, where is the {N + l)-st column of DC 

We show in Proposition 19. 1 1 of the Appendix that Di: 7 v = — JDj.jyJ, where J is 
the exchange matrix with ones on its counterdiagonal and zeros elsewhere. It follows 
that = -J(DLvr)-' J. Consequently, the elements in are the negative 

of the elements in (D|,^)“^, but rearranged. As a result, (D|.^)“^ also possesses 
properties (PI) and (P2), and the analysis of the discrete costate closely parallels the 
analysis of the state. The main difference is that the costate equation contains the 
additional A^f+i term along with the additional equation (16.6() . By ()6.6() and the 
previously established bound ||X||oo < c||y||oo, it follows that 

IIAtv+iIIoo < c||y||oo, (6.7) 

where c is independent of N. Since Dtl = 0, we deduce that (dJ.^)“^Djv+i = —1. 
It follows that 

(Dt)-iD]Y^^ = [(Dj,jv)”^Oln][D]v+iOl»] = -l®In. 

Exploiting this identity, the analogue of (lO) is 

A = (I + (Dt)-iAT)-i[(l 0 I„)Avr+i + (Dt)-i(y3 - QX - SU)]. 

Hence, we have 

||A||oo < 2||(1 0 I„)AAr+i + (D^)-'(y3 - QX - SU)||oo. 

Moreover, ||(1 0 I„)A7v+i||oo_< c||y||oo by (jUTj) and ||(Dt)“V3||oo < 2||y3j|oo- The 
terms || (Dt)“^QX||oo and || (Dt)~^SU)||oo are handled exactly as the term ||D~^BU||oo 
was handled in the state equation dSH). We again conclude that ||A||oo < c||y||oo 
where c is independent of N. 
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Finally, let us examine the fifth component of the linear system VT*[X,U, A] = 
y. These equations can be written 

+ RjUi + = y5i, l<i<N. 

By (A2) the smallest eigenvalue of R,; is greater than a > 0. Consequently, the bounds 
||X||oo < c||y||oo and || A||oo < c||y||oo hnply the existence of a constant c, independent 
of N, such that ||U||oo < c||y||oo. In summary, we have the following result: 

Lemma 6.1. If (A2)-(A3) and (P1)-(P2) hold, then there exists a constant c, 
independent of N, such that the solution o/V7”*[X, U, A] = y satisfies 

||X||oo + ||U||oo + ||A|U<c||y||oo. 


Let us now prove Theorem 11.11 using Proposition 12.11 By Lemma 16.11 p = 
llvr(x*,u*,A*)-i||oo is bounded uniformly in N. Choose e small enough that ep < 
1. When we compute the difference Vr(X, U, A) - Vr(X*,U*,A*) for (X, U, A) 
near (X*,U*, A*) in the oo-norm, the D and constant terms cancel, and we are 
left with terms involving the difference of derivatives of f or C up to second order at 
nearby points. By assumption, these second derivative are uniformly continuous on 
the closure of fl and on a ball around x*(l). Hence, for r sufficiently small, we have 

||Vr(X,U,A)-vr(X*,U*,A*)||oo <e 

whenever 

max{||X - X*||oo, ||U - U*|U, ||A - A*|U} < r. (6.8) 

By Lemma 13.11 it follows that ||T(X*,U*, A*)|| < (1 — pe)r/p for all N sufficiently 
large. Hence, by ProDOsition l2.ll there exists a solution to T(X, U, A) = 0 satisfying 
(16.811 . Moreover, by (I2J1) and the estimate (11.1011 holds. To complete the proof, 
we need to show that (X,U) is a local minimizer for (12.111 . After replacing the KKT 
multipliers by the transformed quantities given by (1^ . the Hessian of the Lagrangian 
is the following block diagonal matrix: 

diag{cuiV^,,_„)R(Xi,Ui,Ai), ... , V2C'(XAr+i)} 

where H is the Hamiltonian. In computing the Hessian, we assume that the X and 
U variables are arranged in the following order: Xi, Ui, X 2 , U 2 , ..., Xjv, Uat, 
Xjv+i. By (A2) the Hessian is positive definite when evaluated at (X*,U*,A*). 
By continuity of the second derivative of C and f and by the convergence result 
(11.1011 . we conclude that the Hessian of the Lagrangian, evaluated at the solution of 
T(X, U, A) = 0 satisfying (j6.811 . is positive definite for N sufficiently large. Hence, 
by the second-order sufficient optimality condition [20j Thm. 12.6], (X,U) is a strict 
local minimizer of dm. This completes the proof of Theorem 11.11 

7. Numerical illustration. Although the assumptions (A1)-(A3) are sufficient 
for exponential convergence, the following example indicates that these assumptions 
are conservative. Let us consider the unconstrained control problem 

inin{—a:(2) : x{t) = |(—a;(t) + x{t)u{t) — u(t)‘^), x(0) = l} . (7.1) 
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Fig. 7.1. The base 10 logarithm of the error in the sup-norm as a function of the number of 
collocation points. 


The optimal solution and associated costate are 

X* {t) = 4:/a{t), a(t) = 1 + 3 exp(2.5t), 
u*{t)=x*{t)/2, 

\*{t) = exp(—2.5t)/[exp(—5) + 9exp(5) + 6]. 

Figure rm plots the logarithm of the sup-norm error in the state, control, and costate 
as a function of the number of collocation points. Since these plots are nearly linear, 
the error behaves like clO““^ where a ~ 0.6 for either the state or the control and 
a ~ 0.8 for the costate. In Theorem ll.il the dependence of the error on N is somewhat 
complex due to the connection between m and N. As we increase TV, we can also 
increase m when the solution is infinitely differentiable, however, the norm of the 
derivatives also enters into the error bound as in (lO) . Nonetheless, in cases where 
the solution derivatives can be bounded by c"* for some constant c, it is possible to 
deduce an exponential decay rate for the error as observed in m Sect. 2]. Note that 
the example problem dm does not satisfy (A2) since V^C = 0 , which is not positive 
definite. Nonetheless, the pointwise error decays exponentially fast. 

8. Conclusions. A Gauss collocation scheme is analyzed for an unconstrained 
control problem. For a smooth solution whose Hamiltonian satisfies a strong convexity 
assumption, we show that the discrete problem has a local minimizer in a neighbor¬ 
hood of the continuous solution, and as the number of collocation points increases, the 
distance in the sup-norm between the discrete solution and the continuous solution is 
0{N'^~'^) when the continuous solution has rj -\- 1 continuous derivatives, rj > 3, and 
the number of collocation points N is sufficiently large. A numerical example is given 
which exhibits an exponential convergence rate. 

9. Appendix. In (12. 4|) we define a new matrix in terms of the differentiation 
matrix D. The following proposition shows that the elements of are the negative 
and a rearrangement of the elements of D. 
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Proposition 9.1. The entries of the matrices D and satisfy 
Dij = —1 < * < ^ < j < N. 

In other words, Di-jv = —where J is the exchange matrix with ones on its 
counterdiagonal and zeros elsewhere. 

Proof. By (11.91) the elements of D can be expressed in terms of the derivatives of 
a set of Lagrange basis functions evaluated at the collocation points: 

A, = L,{n) where L, € L,{n) = | J q < fc < IV, M J- 

In (11.91) we give an explicit formula for the Lagrange basis functions, while here we 
express the basis function in terms of polynomials Lj that equal one at tj and vanish 
at Tfc where Q<k<N,k^j. These N + 1 conditions uniquely define Lj € Vn. 
Similarly, by [TU Thm. 1], the entries of are given by 

4 =M,(r,) where M, (r,) = | J ti2i<N + l,kj^j. 

Observe that MM+i-j{t) = Lj{—t) due the symmetry of the quadrature points around 
t = 0: 

(a) Since —tn+i-j = Tj, we have Lj{—TN+i-j) = LjiTj) = 1- 

(b) Since tn+i = 1 and tq = —1, we have Lj(—r^v+i) = Lj^tq) = 0. 

(c) Since —Ti = tn+i-i, we have Lj{—Ti) = Lj{TN+i-i) = 0 if i 7 ^ iV + 1 — j. 
Since M^+i-jif) is equal to Lj{—t) at + 1 distinct points, the two polynomials are 
equal everywhere. Replacing MN+i-j[t) by Lj{—t), we have 

-^jv+i-i.Ar+i-j = —Lj[—T]si+i-i) = —Lj{Ti) = —Dij. 


□ 

Tables I^TT] and 1^1^ illustrate properties (PI) and (P2) for the differentiation matrix 
D. In Table [TT] we observe that ||D(l]y||oo monotonically approaches the upper limit 
2. More precisely, it is found that ||Djl]y||oo = 1 +tn, where the final collocation 
point tn approaches one as N tends to inhnity. In Table 19.21 we show the maximum 
2-norm of the rows of It is found that the maximum is attained by 

the last row of and the maximum monotonically approaches v^- 


N 

25 

50 

75 

100 

125 

150 

norm 

1.995557 

1.998866 

1.999494 

1.999714 

1.999816 

1.999872 

N 

175 

200 

225 

250 

275 

300 

norm 

1.999906 

1.999928 

1.999943 

1.999954 

1.999962 

1.999968 


Table 9.1 
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N 

25 

50 

75 

100 

125 

150 

norm 

1.412201 

1.413703 

1.413985 

1.414085 

1.414131 

1.414156 

N 

175 

200 

225 

250 

275 

300 

norm 

1.414171 

1.414181 

1.414188 

1.414193 

1.414196 

1.414199 


Table 9.2 

Maximum Euclidean norm for the rows of [W^/^Di;iv]~^ 
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